;; Non-robust computational geometry functions
;; Copyright (C) 2021  Rocks Mazama
;;
;; This program is free software: you can redistribute it and/or modify
;; it under the terms of the GNU General Public License as published by
;; the Free Software Foundation, either version 3 of the License, or
;; (at your option) any later version.
;;
;; This program is distributed in the hope that it will be useful,
;; but WITHOUT ANY WARRANTY; without even the implied warranty of
;; MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
;; GNU General Public License for more details.
;;
;; You should have received a copy of the GNU General Public License
;; along with this program.  If not, see <https://www.gnu.org/licenses/>.

(ns mazama.compgeom.segment-intersection
  (:require [mazama.compgeom.core :refer [interval-intersect?
                                          perp-dot-product
                                          point-orient2d]]))

;;;; NOTE
;;;;
;;;; 1. This software may fail if floating point values are used.
;;;;    This is because == is used to compare numbers and no special effort
;;;;    is made to determine if two floating point values are nearly equal.
;;;;
;;;; 2. This software assumes that:
;;;;    the x-axis increases to the right;
;;;;    the y-axis increases DOWNWARD.

(defn- ratio-unit-inverval?
  "Return true if and only if 0 <= n/d <= 1.
  Both n and d may be either negative, zero, or positive."
  [n d]
  ;; We accomplish this without actually doing the division
  (or (and (pos? d) (<= 0 n d))
      (and (neg? d) (<= d n 0))))

(defn- colinear-overlap-result []
  {:code :colinear-overlap})

(defn- proper-intersection-result [point]
  {:code :proper
   :point point})

(defn- vertex-intersection-result [point]
  {:code :vertex
   :point point})

(defn- segment-parallel
  "Handle the segment intersection case when the two segments are parallel.
  *Note*: It is assumed that the bounding boxes of the segments are known
  to intersect."
  [x1 y1 x2 y2 x3 y3 x4 y4]
  ;; Since the segments are parallel they must be colinear to intersect
  (when (and (= 0 (point-orient2d x1 y1 x2 y2 x3 y3))
             (= 0 (point-orient2d x2 y2 x3 y3 x4 y4)))
    ;; The segments are colinear and, since we are assuming the bounding
    ;; boxes intersect, the segments either intersect at a single endpoint
    ;; or overlap.
    (let [z (cond
              (and (== x1 x3)
                   (== y1 y3)
                   (<= (* (- x4 x1) (- x2 x1)) 0)
                   (<= (* (- y4 y1) (- y2 y1)) 0))
              [x1 y1]
              (and (== x1 x4)
                   (== y1 y4)
                   (<= (* (- x3 x1) (- x2 x1)) 0)
                   (<= (* (- y3 y1) (- y2 y1)) 0))
              [x1 y1]
              (and (== x2 x3)
                   (== y2 y3)
                   (<= (* (- x1 x2) (- x4 x2)) 0)
                   (<= (* (- y1 y2) (- y4 y2)) 0))
              [x2 y2]
              (and (== x2 x4)
                   (== y2 y4)
                   (<= (* (- x1 x2) (- x3 x2)) 0)
                   (<= (* (- y1 y2) (- y3 y2)) 0))
              [x2 y2])]
      (if z
        (vertex-intersection-result z)
        ;; The segments intersect at a non-endpoint. If one of the segments
        ;; is degenerate, it is considered a proper intersection (since the
        ;; the intersection is a single point), otherwise it is an overlap.
        (let [z (cond
                  (and (== x1 x2)
                       (== y1 y2))
                  [x1 y1]
                  (and (== x3 x4)
                       (== y3 y4))
                  [x3 y3])]
          (if z
            (proper-intersection-result z)
            (colinear-overlap-result)))))))

(defn- segment-bb-intersect?
  "Do the bounding boxes of the two line segments ab and cd intersect?"
  [ax ay bx by cx cy dx dy]
  (and (interval-intersect? (min ax bx) (max ax bx)
                            (min cx dx) (max cx dx))
       (interval-intersect? (min ay by) (max ay by)
                            (min cy dy) (max cy dy))))

;;; References
;;; https://en.wikipedia.org/wiki/Line%E2%80%93line_intersection
;;; Computational Geometry in C, 2nd edition, section 7.2, code segseg.c
;;; Faster Line Segment Intersection, Franklin Antonio, Graphic Gems III
;;; The Pleasures of "Perp Dot" Products, Example 2, F.S. Hill, Jr., Graphics Gems IV
(defn segment-intersection
  "Compute the intersection of two line segments, ab and cd, defined by
  the points a, b, c, and d.
  One of the following values is returned:
  nil: If the segments do not intersect.
  {:code :vertex, :point [x y]}: If the segments intersect only at an
  endpoint [x y].
  {:code :proper, :point [x y]}: If the segments intersect at a point [x y]
  other than an endpoint.
  {:code :colinear-overlap}: If the segments are colinear and intersect at
  more than just a single point."
  ([[x1 y1] [x2 y2] [x3 y3] [x4 y4]]
   (segment-intersection x1 y1 x2 y2 x3 y3 x4 y4))
  ([x1 y1 x2 y2 x3 y3 x4 y4]
   ;; Ensure the bounding boxes of the two line segments intersect
   (when (segment-bb-intersect? x1 y1 x2 y2 x3 y3 x4 y4)
     (let [AX (- x2 x1)
           AY (- y2 y1)
           BX (- x4 x3)
           BY (- y4 y3)
           D (perp-dot-product BX BY AX AY)]
       (if (zero? D)
         ;; The segments are parallel
         (segment-parallel x1 y1 x2 y2 x3 y3 x4 y4)
         ;; The segments are not parallel
         (let [CX (- x3 x1)
               CY (- y3 y1)
               T (perp-dot-product BX BY CX CY)]
           (when (ratio-unit-inverval? T D)
             (let [U (perp-dot-product AX AY CX CY)]
               (when (ratio-unit-inverval? U D)
                 (let [alpha (/ T D)
                       x (+ x1 (* alpha AX))
                       y (+ y1 (* alpha AY))]
                   (if (and (or (zero? T) (== T D))
                            (or (zero? U) (== U D)))
                     (vertex-intersection-result [x y])
                     (proper-intersection-result [x y]))))))))))))
